/*   Copyright © 2009 Mark Gil Torres
 *
 *   This file is part of mtoolkit.
 *
 *   mtoolkit is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   any later version.
 *
 *   mtoolkit is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with mtoolkit.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef MATRIX_H
#define MATRIX_H

#include <iostream>
#include <valarray>

namespace mt {

  template <typename T> class Matrix;
  template <typename T> Matrix<T> const operator+(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
  template <typename T> Matrix<T> const operator-(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
  template <typename T> Matrix<T> const operator*(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
  template <typename T> Matrix<T> const operator*(T const& inValue1, Matrix<T> const& inValue2);
  template <typename T> Matrix<T> const transpose(Matrix<T> const& inValue);
  template <typename T> Matrix<T> const IdentityMatrix(T const& inValue);

  template <typename T = double>
  class Matrix {
    public:
      Matrix(void);
      Matrix(int const& inValue1, int const& inValue2, std::valarray<T> const& inValue3);
//assignment operator
      Matrix<T>& operator=(Matrix<T> const& inValue);
//index operator
      T& operator[](size_t inValue);
      T& operator()(int const& inValue1, int const& inValue2);
//getter
      int& getRows(void);
      int& getColumns(void);
      std::valarray<T>& getMValue(void);
//setter
      void setMValue(int const& inValue1, int const& inValue2, std::valarray<T> const& inValue3);
      void setRows(int const& inValue);
      void setColumns(int const& inValue);
//matrix operators
      friend Matrix<T> const operator+ <>(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
      friend Matrix<T> const operator- <>(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
      friend Matrix<T> const operator* <>(Matrix<T> const& inValue1, Matrix<T> const& inValue2);
//scalar multiplication
      friend Matrix<T> const operator* <>(T const& inValue1, Matrix<T> const& inValue2);
//transpose operation
      friend Matrix<T> const transpose<>(Matrix<T> const& inValue);
//ready-made matrix
      friend Matrix<T> const IdentityMatrix<>(T const& inValue);
//info member functions
      int const size(void);

    private:
      int rows;
      int columns;
      std::valarray<T> mValue;
  };

  template <typename T>
  Matrix<T>::Matrix(void) {
    rows = 0;
    columns = 0;
  }

  template <typename T>
  Matrix<T>::Matrix(int const& inValue1, int const& inValue2, std::valarray<T> const& inValue3) {
    rows = inValue1;
    columns = inValue2;
    mValue.resize(rows*columns);
    mValue = inValue3;
  }
//assignment operator implementation
  template <typename T>
  Matrix<T>& Matrix<T>::operator=(Matrix<T> const& inValue) {
    rows = inValue.rows;
    columns = inValue.columns;
    mValue.resize(rows*columns);
    mValue = inValue.mValue;
    return *this;
  }

//index operator implementation
  template <typename T>
  T& Matrix<T>::operator[](size_t inValue) {
    return mValue[inValue];
  }

  template <typename T>
  T& Matrix<T>::operator()(int const& inValue1, int const& inValue2) {
    int temp;
    temp = (inValue1 * rows) + inValue2;
    return mValue[temp];
  }
//getter implementation
  template <typename T>
  int& Matrix<T>::getRows(void) {
    return rows;
  }

  template <typename T>
  int& Matrix<T>::getColumns(void) {
    return columns;
  }

  template <typename T>
  std::valarray<T>& Matrix<T>::getMValue(void) {
    return mValue;
  }

//setter implementation
  template <typename T>
  void Matrix<T>::setMValue(int const& inValue1, int const& inValue2, std::valarray<T> const& inValue3) {
    rows = inValue1;
    columns = inValue2;
    mValue.resize(rows*columns);
    mValue = inValue3;
  }

  template <typename T>
  void Matrix<T>::setRows(int const& inValue) {
    rows = inValue;
  }

  template <typename T>
  void Matrix<T>::setColumns(int const& inValue) {
    columns = inValue;
  }

//matrix operators implementation
  template <typename T>
  Matrix<T> const operator+(Matrix<T> const& inValue1, Matrix<T> const& inValue2) {
    Matrix<T> temp;
    int tmpInValue1R = inValue1.rows;
    int tmpInValue1C = inValue1.columns;
    int tmpInValue2R = inValue2.rows;
    int tmpInValue2C = inValue2.columns;
    if (tmpInValue1R == tmpInValue2R && tmpInValue1C == tmpInValue2C) {
      temp.setMValue(tmpInValue1R, tmpInValue1C, (inValue1.mValue + inValue2.mValue));
    }
    else {
      std::cout << "ERROR : The given Matrices are NOT conformable for addition.\n";
    }
    return temp;
  }

  template <typename T>
  Matrix<T> const operator-(Matrix<T> const& inValue1, Matrix<T> const& inValue2) {
    Matrix<T> temp;
    int tmpInValue1R = inValue1.rows;
    int tmpInValue1C = inValue1.columns;
    int tmpInValue2R = inValue2.rows;
    int tmpInValue2C = inValue2.columns;
    if (tmpInValue1R == tmpInValue2R && tmpInValue1C == tmpInValue2C) {
      temp.setMValue(tmpInValue1R, tmpInValue1C, (inValue1.mValue - inValue2.mValue));
    }
    else {
      std::cout << "ERROR : The given Matrices are NOT conformable for subtraction.\n";
    }
    return temp;
  }

  template <typename T>
  Matrix<T> const operator*(Matrix<T> const& inValue1, Matrix<T> const& inValue2) {
    Matrix<T> temp;
    int tmpInValue1R = inValue1.rows;
    int tmpInValue1C = inValue1.columns;
    int tmpInValue2R = inValue2.rows;
    int tmpInValue2C = inValue2.columns;
    if (tmpInValue1C == tmpInValue2R) {
      int tempSize = tmpInValue1R * tmpInValue2C;
      temp.setRows(tmpInValue1R);
      temp.setColumns(tmpInValue2C);
      temp.mValue.resize(tempSize);
      for (int i = 0; i < tempSize; i++) {
        int startIterJ;
        if (tmpInValue1R == tmpInValue2R) startIterJ = (i/tmpInValue1C) * tmpInValue1C;
        else if (tmpInValue1R > tmpInValue2R) startIterJ = (i/(tmpInValue1C+1)) * tmpInValue1C;
        else startIterJ = (i/tmpInValue1C) * tmpInValue1C;
        int endIterJ = startIterJ + tmpInValue1C;
        int c = 0;
        for (int j = startIterJ; j < endIterJ; j++) {
          int cIndex = j;
          int tIndex = ((j - startIterJ) * tmpInValue2C) + (i % tmpInValue2C);
          c += inValue1.mValue[cIndex] * inValue2.mValue[tIndex];
        }
        temp[i] = c;
      }
    }
    else {
      std::cout << "ERROR : The given Matrices are NOT conformable for multiplication.\n";
    }
    return temp;
  }

//scalar multiplication operation
  template <typename T>
  Matrix<T> const operator*(T const& inValue1, Matrix<T> const& inValue2) {
    Matrix<T> temp;
    int tmpInValue2R = inValue2.rows;
    int tmpInValue2C = inValue2.columns;
    temp.setMValue(tmpInValue2R, tmpInValue2C, (inValue1 * inValue2.mValue));
    return temp;
  }

//transpose operation
  template <typename T>
  Matrix<T> const transpose(Matrix<T> const& inValue) {
    Matrix<T> temp;
    int tmpInValueR = inValue.rows;
    int tmpInValueC = inValue.columns;
    temp.setRows(tmpInValueC);
    temp.setColumns(tmpInValueR);
    temp.mValue.resize(tmpInValueC * tmpInValueR);
    for (int i = 0; i < tmpInValueC; i++) {
      for (int j = 0; j < tmpInValueR; j++) {
        int cIndex = (tmpInValueR * i) + j;
        int tIndex = (tmpInValueC * j) + i;
        temp[cIndex] = inValue.mValue[tIndex];
      }
    }
    return temp;
  }

//ready-made matrix
  template <typename T>
  Matrix<T> const IdentityMatrix(T const& inValue) {
    Matrix<T> temp;
    temp.setRows(inValue);
    temp.setColumns(inValue);
    int tempSize = inValue * inValue;
    temp.mValue.resize(tempSize);
    for (int i = 0; i < tempSize; i++) temp[i] = 0;
    for (int j = 0; j < tempSize; j+=(inValue+1)) temp[j] = 1;
    return temp;
  }

//info member functions
  template <typename T>
  int const Matrix<T>::size(void) {
    return (rows * columns);
  }

}

#endif
